Instabilities leading to vortex lattice formation in rotating Bose-Einstein condensates 
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We present a comprehensive theoretical study of vortex lattice formation in atomic Bose-Einstein 
condensates confined by a rotating elliptical trap. We consider rotating solutions of the classical 
hydrodynamic equations, their response to perturbations, as well as time-dependent simulations. 
We discriminate three distinct, experimentally testable, regimes of instability: ripple, interbranch, 
and catastrophic. Under symmetry-breaking perturbations these instabilities lead to lattice for- 
mation even at zero temperature. While our results are consistent with previous theoretical and 
experimental results, they shed new light on lattice formation. 
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Vortex lattices are a striking manifestation of super- 
fluidity in rotating systems. In recent years such states 
have been generated in atomic Bose-Einstein conden- 
sates (BECs) via rotation in an elliptical harmonic trap 
[H Qi y| I in analogy to the classic rotating bucket exper- 
iment in superfluid Hehum jj]. Vortex lattices form for 
trap rotation frequencies fi in the region of SI ~ Q.7uj±, 
where uj± is the mean harmonic trap frequency in the 
rotating plane. Although a vortex lattice is thermody- 
namically favourable for much lower rotation frequen- 
cies, the instabilities necessary to seed vortex lattice for- 
mation are only present in the region of O ~ 0.7uj±. 
Such instabilities have been predicted by hydrodynamic 
studies of condensate solutions in the rotating frame Q 
and the dynamical perturbations of these states , and 
have also been analysed in dipolar BECs [3- Although 
numerical simulations of the Gross-Pitaevskii equation 
have also observed vortex lattice formation in this region 
mmpjIlllElIll, the results are mixed: for Q < 0.7uj±, 
Refs. la S [l3 require thermal effects to enable vortex 
nucleation while Ref. [ll| does not; for fJ^O.Twj,, Refs. 
[3, llil U2, nil observe a shape instability before nucleat- 
ing vortices even in the absence of thermal effects. This 
is consistent with experiments 0, [iJI which indicate that 
lattice formation is temperature-independent. Breaking 
the rotational symmetry of such simulations has been 
shown to be crucial to enable realistic nucleation of vor- 
tices o. 

In this paper we present a thorough investigation of 
the instabilities leading to vortex lattice formation in 
elliptical traps in the regime f2 < u!±. We primarily 
consider the cases when the trap ellipticity or rotation 
frequency is introduced adiabatically, for which it is ap- 
propriate to transform to the rotating frame. In order to 
probe the condensate instabilities we consider the rotat- 
ing frame condensate solutions within the classical hy- 
drodynamic regime and the response of these states to 
dynamical perturbations. Furthermore we perform time- 
dependent (TD) condensate simulations to test our re- 
sults, and probe the fate of the unstable condensate. 

In the experiments [l|, y| the harmonic trap confining 
the BEG is transformed radially into an ellipse and ro- 



tated about the z-axis. In the limit of zero temperature, 
the condensate can be approximated by a mean-field 
'wavefunction' ip, which can be expressed as '0 = yfpe'^'^, 
where p and (f> are the condensate density and phase. In 
the frame rota ting at f2, the density p and fluid velocity 
V ~ {h/m)W(j) [l5j satisfy the hydrodynamic equations. 
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Here m is the atomic mass and g ~ AttTi a/m is the inter- 
action coefficient, where a is the s-wave scattering length 
(in this work g > 0). The harmonic trap is defined by 
the frequencies uJx , ujy and lOz ■ When the trap is elliptical 
in the radial plane we can write the x, y frequencies as 
ujx = -\/l — e Lo± and LOy — \J\ -|- e uj^, where e charac- 
terises the trap ellipticity, and lo\ = (w^ -I- a;y)/2. 

Following the approach by Recati et al. Q , we employ 
the Thomas- Fermi (TF) approximation by neglecting the 
^'^ \fp I \fp) term in Eq. (01. Furthermore we assume an 
irrotational quadrupolar velocity field |l6| , 



V = aSIixy), 



(3) 



where a is the velocity field amplitude. Substituting this 
into Eq. j^jl and setting dp/dt = dv/dt = 0, we obtain 
stationary density solutions of the form, 
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in the region where the chemical potential p, > m{u)'^x'^ + 
tjyj/^ -I- a;^z^)/2, and po = elsewhere. The effect of the 
rotation is to give effective trap frequencies, which satisfy, 
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FIG. 1: (a) Velocity field amplitude a of the stationary so- 
lutions of Eq. (0 as a function of rotation frequency fi for 
e = (solid line), 0.02 (dotted line) and 0.1 (dashed line). Re- 
gions of dynamical instability for e = 0.02 and 0.1 are shown 
(circles), (b) Phase diagram of e versus f2. Plotted are the 
bifurcation point f2bif(e) (solid line) from Eq. l|7|l, the on- 
set of dynamical instability f2ins(e) from Eq. (|HJ (dashed line), 
and experimental data of Hodby ei al. [3 (circles). Further- 
more, time-dependent simulations of Eqs. Q and Q show 
the critical ellipticities beyond which the condensate deviates 
from an elliptical shape e^J^^ (crosses) and beyond which lat- 
tice formation ultimately occurs ej.""* (points with error bars) . 
The bifurcation point for e = (dotted line) and crossing fre- 
quency fix (dot-dashed line) at which r2bif(e) = nina(e) are 
indicated. We assume a condensate with ^ = IQfiuj^. 



where a determines the cllipticity of the BEC density Q . 
Plugging Eq. gj into Eq. iQJ we obtain , 



Zjl{a + n)+ul{a-VL) =0, 



(7) 



which specifies the stationary condensate solutions in the 
frame rotating at fJ (in the TF regime). 

Figure 1(a) shows the stationary solutions in ^-a space 
for various values of trap ellipticity e. In the limit of 
e = (solid line) a non-rotating (a = 0) solution occurs 
for all ri, with two additional solutions bifurcating from 
the a = axis at fij^jf = u;_l/-\/2. For finite e (dotted 
and dashed lines), the a = solution disappears and the 
plot consists of two distinct branches. The upper branch 
(positive a) is single- valued and exists for all il, while 
the lower branch (negative a) is double-valued and ex- 
ists only when f2 is greater than the bifurcation frequency 
51bif(e)- As e is increased from zero, the branches move 
away from the a = axis, as can be seen in Fig. ^a). 
Furthermore the bifurcation point of the lower branch 
^bif(e) shifts to higher rotation frequencies, as shown in 
Fig- db) (solid line). The branch diagrams, which have 
been probed experimentally [3 , are key to understanding 



the response of the system to the adiabatic introduction 
of trap ellipticity e or rotation frequency fi. Before any 
rotation/ellipticity is applied the BEC has a = 0. When 
Q, is increased adiabatically (for fixed e), the BEC fol- 
lows the upper branch, with increasing a and ellipticity 
in the density profile. When e is introduced adiabatically 
(for fixed 17), the BEC can follow two routes, depend- 
ing on the value of f2 relative to the bifurcation point 
fij^if (dotted line in Fig.IHb)). For Vl < l^gjf, the lower 
branch is nonexistent and the BEC follows the upper 
branch to increasing values of a. For Vt > il^jf , the lower 
branch moves from a = to negative a, and the BEC 
follows this route. However, as e is increased, the edge 
of the lower branch fibif (e) shifts to higher frequencies, 
and when ribif(e) > ^ the lower branch no longer exists. 
In this manner, the evolution of the branches can induce 
instability, and has been linked to lattice formation [2|- 

The stationary solutions of Eq. {TJ are not necessarily 
stable solutions. To investigate their stability we follow 
the approach of Sinha and Castin |y| by considering small 
perturbations Sp and S(j) of stationary solutions of den- 
sity po smd phase (po- Taking variational derivatives of 
Eqs. ^ and (j2Jl we obtain the time evolution equations, 
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where Vc = v — Jl x r is the velocity field in the rotat- 
ing frame. The eigcnfunctions of Eq.JS]) grow in time as 
e^*, where A is the corresponding eigenvalue. Solutions 
are unstable when there exist one or more eigcnfunctions 
for which A contains a positive real part. Such unstable 
solutions are thought to seed vortex lattice formation |g] . 

Using a polynomial ansatz for the perturbations, this 
method predicts a region of dynamical instability when 
fi exceeds a critical value fiins(e) 13, as indicated in 
Fig- da) (circles). The unstable modes have lengthscales 
of the order of the condensate size, much greater than 
the healing length ^ = h/y^mnog which characterises the 
vortex size. In the limit of e = 0, f2ins ~ O.TSw^. As e 
is increased, i7ins(e) is reduced (dashed line in Fig.QJb)). 
Note that the outlying point in Fig. d^) for e = 0.1 
(dashed line) at 17 w 0.61uj± is not considered to be in the 
region of instability due to its narrow width and compar- 
atively small eigenvalues [I3I- Note that, for the param- 
eter space of interest, on the lower branch the solution 
closest to the a = axis is never dynamically unstable. 
A key rotation frequency in our work is fix, which is the 
crossing point of f2bif(e) and f7iiis(e), and has the value 
fix w 0.765uj±_ (dot-dashed hne in Fig.db)). 

Based on the stationary solutions of Eq. {TJ and the dy- 
namical instability of Eq. (O we can predict the stability 
of a BEC following the adiabatic introduction of fl or e. 
However, to determine how the instability manifests itself 
and whether it leads to lattice formation we perform TD 
simulations of Eqs. d and id in the laboratory frame. 
This is equivalent to solving the Gross-Pitaevskii equa- 
tion [T3| ■ Following previous approaches [^ [13, Ull H^ , 
and noting that the solutions of Eqs. d and © are 



independent of z [13 , we consider a 2D system. The ini- 
tial state is found by imaginary-time propagation. Then, 
in real time, either e is ramped up linearly at a rate of 
de/dt = 10~'*(jjj_ (for fixed f2), or f2 is ramped up linearly 



at a rate dfl/dt = 10 



(for fixed e). We monitor the 



evolution of a by fitting the velocity field in a central 
region [3 x 3]fim to the form of Eq. ||21l. 

In such 'idealised' simulations the trap and BEC can 
maintain a two-fold rotational symmetry to unrealistic 
levels. This strongly inhibits vortex nuclcation since 
the vortices must enter in pairs rather than individu- 
ally. Symmetry-breaking, rather than thermal effects, 
has been shown to the crucial in simulating lattice for- 
mation [IjI. Indeed, Eq. ((HJ predicts that only odd modes 
of the system are dynamically unstable. Previous studies 
|8|, lid ] have broken this symmetry through the introduc- 
tion of thermal 'noise'. To overcome this problem we 
shift the trap by half a healing length before running in 
real-time, thereby allowing excitation of odd modes |l7l| . 

We first consider the adiabatic increase of trap ellip- 
ticity e for fixed CI. We discriminate three cases of in- 
stability, which each occur in distinct frequency regimes, 
as indicated in Fig.njb): (I) ripple instability, (II) inter- 
branch instability, and (III) catastrophic instability. Wc 
will discuss each in turn. 

/ Ripple instability, fl < uj±/\/2. The case for O = 
0.7iu± is shown in Fig. [3(a). The velocity field ampli- 
tude a (dots) follows the upper branch of the rotating 
solutions (solid line), for which the condensate axes ro- 
tate in phase with the trap axes, as noted experimentally 
0,l3- The upper branch is always a solution, and, as e 
is increased, the condensate moves along the branch to 
higher a. However, when the ellipticity exceeds a critical 
value ecr (corresponding to when Q > iljnst) the solution 
becomes dynamically unstable, according to Eq. ©. For 
the example in Fig. [3(a), ecr ~ 0.09 (dashed line). Sub- 
sequently a (dots) deviates from the rotating solutions 
of Eq. O (solid lines), consistent with the onset of dy- 
namical instability. For low e, a (dots) features small 
oscillations due to the centre-of-mass motion caused by 
the initial offset of the condensate. 

At a critical ellipticity e'^^^ low density ripples form on 
the condensate edge (Fig. [3(b)), each on the scale of the 
healing length and featuring a phase singularity ('ghost' 
vortices |l,|ajll3)- These ripples grow in amplitude as 
e is increased. If e becomes fixed when the ripples have 
very low amplitude they do not grow over the timescales 
considered. However, once e exceeds a second critical 
value ej."*' (corresponding to when the ripples have am- 
plitude of order of 10%rio), the dynanncal instability of 
Eq. (O is triggered by the ripples. This instability gener- 
ates largescale shape osciUations (Fig. [3(c)), disrupting 
the condensate, and enabling 'ghost' vortices to nucle- 
ate into the condensate, which slowly crystallise into a 
lattice (Fig. [3(d)) HI. Once ej,'f is reached, lattice for- 
mation occurs independently of whether e becomes fixed 
or continuously increased. Since the ripples that trigger 
the dynamical instability originate in the non-TF tails of 
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FIG. 2: Dynamics under a continuous increase of e (at a rate 
of de/dt = lO^^t^i) for I. Ripple instability (r2/a;j_ = 0.7), 11. 
Interbranch instability {Q/uj± — 0.75), and III. Catastrophic 
instability {fl/uj±_ = 0.8). (a) Velocity field amplitude a ver- 
sus e according to Eq. Q (solid lines) and TD simulations of 
Eqs. Q and l[3 (dots). In the latter case, a is derived by 
fitting the velocity field to Eq. ^. To the right of the dashed 
line this is no longer a valid fit (the standard deviation of 
the fit becomes of the order of a). The regions of dynamical 
instability are indicated (circles). Density snapshots corre- 
sponding to (b) the onset of instability, (c) disrupted state 
seeded by the instability, and (d) a vortex lattice. Once the 
instability point is reached, the dynamics are qualitatively 
similar whether e is continuously increased or becomes fixed. 
Dark/light regions represent high/low density. Each box rep- 
resents a region (12 x 12)iJ.m (the simulation box is much 
greater than this). In 1(b) and 11(b) the density scale is lim- 
ited to 0.1%no to highlight low density features. 



the BEC, they cannot be explained by the TF analytics 
of Eqs. (3) - (8), and a higher-order (non-TF) approach 
would be required to explain their origin. 

Surface ripples have been observed experimentally to 
precede vortex nuclcation at this rotation frequency [3]. 
The gradual growth of the ripples leads to a slow injec- 
tion of energy/angular momentum into the system, as 
has been observed in previous studies within this fre- 
quency regime H, [ll|. Note that according to Eq. © 
the dynamical instability on the upper branch couples to 
odd modes only. If the symmetry is preserved we do not 
expect the instability to develop [3, S HOl ■ 

// Interbranch instability, (jJ±/\/2 < f2 < fix- The 
case for il ~ 0.75uj± is shown in Fig. [31(a). Since 
ri > u!j_/V2, a (dots) initially follows the lower branch 
solutions, where we observe that the BEC and trap axes 



are 7r/2 out of phase. As e is increased a point is reached 
when O < J7bif(£)- Here the lower branch is no longer a 
solution of Eq. ^. For the example shown, this occurs 
for e « 0.02 (dashed line in Fig. [31(a)). Due to the non- 
TF nature of the numerical solutions of Eqs. (QJ and |(5J, 
a does not perfectly fit the branch solutions of Eq. (0) . 

When a (dots) reaches the cusp of the lower branch 
it deviates non-adiabatically, as observed experimentally 
[2|. Since f2 < Six, the upper branch is dynamically sta- 
ble at this ellipticity. The condensate tries to transform 
to the upper branch, but without dissipation it cannot 
relax to this state. Instead a (dots) oscillates between 
positive and negative values, and the condensate den- 
sity undergoes quadrupolar shape oscillations between 
an almost circular and highly elliptical shape. If e be- 
comes fixed close to this critical ellipticity, the quadrupo- 
lar shape oscillations are stable. However, if e is increased 
further, the shape oscillations destabilise, with the con- 
densate shedding low density material at its extrema in a 
spiral pattern (Fig.|3l(b)). This situation is closely anal- 
ogous to when the rotation/ellipticity is suddenly turned 
on, with the fate of the condensate being qualitatively 
similar [i3;[i3- The growth of the ejected material grad- 
ually destabilises the condensate (Fig.|3l(c)), leading to 
vortex nucleation and ultimately a lattice (Fig. 131(c)). 
This is fully consistent with the observations in J2|. 

/// Catastrophic instability, SI > Six- The case for 
O = O.StJ^ is shown in Fig. [3n(a). Again a (dotted 
line) follows the lower branch, which ceases to be a so- 
lution at some critical e (dashed line). However, since 
SI > fix, the upper branch is dynamically unstable at 
this point, and no stable solutions exist. The BEC under- 
goes a quick and catastrophic instability, with a (dots) 
deviating rapidly from the rotating solutions of Eq. l [7| l. 
The BEC density profile becomes strongly contorted into 
a spiral shape (Fig. |3n(b)). The arms of the spiral col- 
lapse inwards and trap phase singularities to form vor- 
tices. Energy and angular momentum are very rapidly 
injected into the BEC (in contrast to the gradua l rip ple 
and interbranch instabilities), as observed in |3, lll| fo^' 
this frequency regime. The BEC becomes highly excited 
and turbulent (Fig. |3n(c)), with structure on length- 
scales less than the healing length. Although we observe 
this state to ultimately settle into a lattice (Fig. |3l(c)), 
one may question the validity of our zero-temperature 
condensate simulations for such a 'heated' state. Indeed, 
for S]>0.78w_L, turbulent states, rather than vortex lat- 
tices, were observed experimentally in []|. 

In the TD simulations we have measured two distinct 
critical ellipticites: ej?°^ (crosses in Fig.QIb)) is when, for 
a continuously increasing e, we observe the density to de- 
viate from a smooth ellipse (on the level of 0.1%no); ej."*'* 
(points with error bars in Fig. 0b)) is when, for e ramped 
up to some final value, instability and lattice formation 
ultimately occur. In regime I, typically e^^'^ < e™*** since 
surface ripples are generated which are stable for a nar- 



row range of e, above which they induce instability and 
lattice formation. In regimes II and III, ej?^^ w ecr*^*: indi- 
cating the relatively sudden onset of this instability once 
the density deviates from a smooth ellipse. 

According to the TD simulations, the region above the 
points in Fig. db) is unstable, leading to lattice forma- 
tion. The prediction of Eq. JHJ (dashed fine) gives rea- 
sonable agreement with the TD results in region I, while 
Eq. (7J) gives excellent agreement in region III. Also plot- 
ted in Fig. nfb) are the experimental results of Hodby 
et al. Q (circles). The TD results give good agreement 
with the experimental data throughout, with the agree- 
ment being particularly good in region HI. 

So far we have considered the adiabatic introduction of 
ellipticity for fixed SI. However, the results in Fig. 0b) 
also apply to the case when O is introduced adiabati- 
cally for fixed e. Here the condensate always follows the 
upper branch and so is only prone to the ripple instabil- 
ity when the upper branch solution becomes dynamically 
unstable. According to Eq. © dynamical instability oc- 
curs when SI > r2inst(e) (dashed hue in Fig.QIb)). In |3], 
for e = 0.025, vortex lattice formation was observed to 
occur when Q^0.75uj±, which agrees well with our TD 
simulations and Fig. 0b). 

Equations and 10) are valid in the TF limit of 
ji/hid^ » 1. For the TD simulations of Eqs. 0) and 
(0); which do not make this assumption, we have em- 
ployed an intermediate value of fj,/?iu!± = 10, but have 
tested other values. For less (more) Thomas- Fermi- Zifce 
condensates, the TD results in Fig. 0b) shift downwards 
(upwards); deviating away from (towards) the TF pre- 
dictions of Eqs. 0) and (jSJ). In particular, as ^/fiLU± is 
increased, the TD values of a in Fig. Il-III(a) (dots) 
become closer to the branch solutions (solid lines). 

We have shown that vortex lattice formation is inher- 
ently a two-dimensional and zero temperature effect. We 
have theoretically mapped out the condensate stability 
as a function of rotation frequency f2 and trap ellipticity 
e, with our interpretation being consistent with previ- 
ous experimental and theoretical results. Specifically, for 
fixed SI and adiabatic introduction of e we find three dis- 
tinct regimes of instability - ripple {Q < ujj_/^/2), inter- 
branch {ljJj_/\/2 < rt < six) and catastrophic (SI > Six)- 
Each instability manifests itself in a characteristic man- 
ner, which could be observed experimentally. Ultimately 
in each case the instability seeds vortex lattice formation. 
This formation process not only relies on the presence of 
an instability but is crucially dependent on the break- 
ing of the two-fold rotational symmetry of the system, as 
inevitably occurs experimentally. 
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